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Abstract 

This paper shows how commuting left and right actions of Lie groups on a manifold may be used to com- 
plement one another in a variational reformulation of optimal control problems as geodesic boundary value 
problems with symmetry. In such problems, the endpoint boundary condition is only specified up to the right 
action of a symmetry group. In this paper we show how to reformulate the problem by introducing extra degrees 
of freedom so that the endpoint condition specifies a single point on the manifold. We prove an equivalence 
theorem to this effect and illustrate it with several examples. In finite-dimensions, we discuss geodesic flows 
on the Lie groups SO(3) and SE(3) under the left and right actions of their respective Lie algebras. In an 
infinite-dimensional example, we discuss optimal large-deformation matching of one closed curve to another 
embedded in the same plane. In the curve-matching example, the manifold Emb(S ,1 ,M 2 ) comprises the space 
of closed curves S 1 embedded in the plane R 2 . The diffeomorphic left action Diff(R 2 ) deforms the curve by a 
smooth invertible time-dependent transformation of the coordinate system in which it is embedded, while leaving 
the parameterisation of the curve invariant. The diffeomorphic right action Diff(S' 1 ) corresponds to a smooth 
invertible reparameterisation of the S domain coordinates of the curve. As we show, this right action unlocks 
an important degree of freedom for geodesically matching the curve shapes using an equivalent fixed boundary 
value problem, without being constrained to match corresponding points along the template and target curves 
at the endpoint in time. 

1 Introduction 

In this paper we are concerned with finding geodesies between points on manifolds. The construction of geodesies 
is useful for studying problems on manifolds since they can describe the relationship between two points. Within 
a coordinate patch on a manifold, any point can described relative to a reference point by specifying a direction 
and a length along the geodesic in that direction. This becomes useful for performing statistics on the coordinate 
patch, for example. In this paper we consider problems in which the endpoint of the trajectory is only fixed up to 
the orbit of a Lie group. In low dimensional cases (and we shall describe some examples of these) it is often easy to 
solve these problems by constructing reduced coordinates which do not change under the action of the Lie group. 
However, in many cases it is difficult to construct such coordinates, especially if the problem is to be discretised and 
solved numerically. In this paper we provide a framework that allows one to work with full unreduced coordinates 
on the manifold, by transforming to an equivalent problem which has the endpoint of the trajectory fixed exactly. 



There are many examples of problems where this framework can be applied, but we are motivated by the prob- 
lem of obtaining diffcomorphisms on M 2 which map one embedded curve T A into another embedded curve T B , 
and which minimise a given metric so that they are geodesies in the diffcomorphism group. The aim is to find a 
characterisation of curve T B with respect to curve T A that is independent of parameterisations of the curves. This 
means that we do not specify a priori the point on T A which gets matched to each specific point on T B , and so the 
minimisation is performed over all parameterisations of the curves. In practise the computation is performed using 
a particular parameterised curve q G Emb(S', R 2 ) (where S is the embedded space, for example, the circle for simple 
closed curves). In computing the equations of motion, a conjugate momentum p q G T* Em^S", R 2 ) is constructed, 
and the flow taking the initial curve T A to the final curve T B can be characterised entirely by the initial conditions 
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p q \t=o for the conjugate momentum. In fact, it turns out that p q \t=o is normal to the curve, so the flow can be 
characterised by a one-dimensional signal. Since T* Emb(S', K 2 ) is a linear space, linear statistics can be computed 
on p q \t=o- For example, this may allow one to test the hypothesis that there is a statistical correlation between be- 
tween the shape of the surface of a biological organ, obtained from a medical scan, and future development of disease. 

To discuss the issues further, we formulate the curve matching problem described above, which may be regarded 
as an optimal control problem in the sense of the problems discussed in [BCHMOO, BCMR98 : 

Definition 1 (Curve matching problem). Let q(s;t) be a one-parameter family of parameterised simple closed 
curves in M 2 , with s G [0, 1] being the curve parameter and t G [0, 1] being the parameter for the family. Let u(x; t) 
be a one-parameter family of vector fields onM 2 . Let n be a diffeomorphism of S 1 . We seek q and u which satisfy 

f 1 1 

mm/ dHv-d* 

subject to the constraints 

d 

[Reconstruction relation] — q(s-.t) = u(q(s;t),t), (1) 

Initial state (Template)] q{s;0) = q A (s), (2) 
[Final state (Target)] q(s;l) = q B (n{s)) 1 (3) 

where \\ ■ \\y is the chosen norm which defines the space of vector fields V. 

The solution of this problem describes a geodesic in the diffeomorphism group which takes the simple closed 
curve T A parameterised by q A to the simple closed curve T B parameterised by q B . We represent the shapes of 
simple closed curves as elements of Emb(5 1 , K 2 )/ Diff(S' 1 ), where Diff(S' 1 ) is the group of diffeomorphisms of S 1 . 
However, we do not want to calculate on this space; instead, we want to calculate on the full space Emb(jf? , M. 2 ) by 
minimising over all reparameterisations n(s) G Diff(S' 1 ). 

There are two general strategics for solving such problems. The first strategy, used for example in |CY01j . 
is to use a gradient method (i.e. a modification of the steepest descent method such as the nonlinear conjugate 
gradient method [She941 and references therein]) to minimise the action integral over paths q(s,t) which satisfy the 
dynamical constraint (this constraint was enforced "softly" via a penalty term in |CY01] L An alternative method, 
referred to in |MMS06| as the "Hamiltonian method", is to introduce Lagrange multipliers p(s,t) which enforce the 
dynamical constraint, and to derive Hamilton's canonical equations for q and p q7 following the general derivation 
described in CH09I for example]. Minimisation over the reparameterisation n, together with a conservation law 
obtained from Noether's theorem, results in the condition that the tangential component of p q vanishes. The aim 
of the Hamiltonian method is to turn an optimisation problem into an algebraic equation given by the time-1 flow 
map of Hamilton's canonical equations. One then solves a shooting problem to find initial conditions for the normal 
component p q which generate solutions to Hamilton's equations that satisfy the boundary condition ([3j). The 
difficulty in solving this problem numerically lies in finding a good numerical discretisation of the target constraint 
condition ([3]). Various functionals have been proposed which vanish when the constraint condition is satisfied. In 
[GTY04) a functional was proposed based on singular densities (measures) , and in [VG05] a functional was proposed 
based on singular vector fields (currents). An alternative spatial discretisation for the current functional based on 
particle- mesh methods was proposed in |Cot08| . There are several difficulties with these functionals: one is that 
after numerical discretisation the functionals do not vanish at the minima, and the boundary condition must be 
replaced by a functional minimising condition. It is also difficult to express the probability distribution of the 
functional given the distribution of measurement errors; this is important for statistical modelling. 

In this paper we consider a transformation of problems of the above type, which results in an alternative for- 
mulation that removes the reparameterisation variable rj from the target constraint, thereby resulting in a standard 
two point boundary value problem on T* Emb(S', R 2 ) (with a constraint on the initial conditions plus an additional 
parameter). This transformation can be applied to a very general class of problems; so we present it in the general 
case of Lie group actions on a manifold. 

The rest of this paper is organised as follows. In Section [2l we formulate the optimal control problem, then 
transform to the geodesic problem with symmetry and prove that the two problems are equivalent. In Section [3] 
we give some examples and discuss the application to matching curves and surfaces. Section 2] is the summary and 
outlook. 
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2 Reparameterised geodesic boundary value problems with symmetry 

In this section we describe a general framework for geodesic boundary value problems with symmetry. We define 
the following Optimal Control Problem. 

Definition 2 (Geodesic boundary value problem with symmetry). Let Q be a manifold, let G be a Lie group 
acting on Q from the left, and let H be a (possibly different) Lie group acting on Q from the right that commutes with 
the left action of G on Q, with corresponding Lie algebras g and \), and corresponding Lie algebra actions X G and 
X H respectively. Furthermore, let A : g — > g* be a positive-definite self-adjoint operator and let (■, -} g : g X g* — > M 
be a nondegenerate pairing which defines an inner product on g. We seek 

• a one parameter family q of points on Q parameterised by t S [0, 1], 

• a one parameter family £ of elements of g for t € [0, 1], and 

• T] £ H, 



which minimise 



subject to the constraints 



[Reconstruction relation] -^q = X G q, (4) 

[Initial state (Template)] q\t=o = q A , (5) 
[Final state (Target)] q\ t=1 = R n q B , (6) 

where q A , q B are chosen points on Q, and R n is the right-action of rj on Q. 

Remark 3. This problem is an optimal control problem in which we seek the shortest path in Q from q A to any 
point q B rj, rj £ H . This means we are seeking the shortest path in Q/H , but are performing the computation on Q. 
In many cases it is much easier to compute on Q, for example when Q is a vector space. We refer to this process 
of solving a problem on Q/H by calculating on Q as "un-reduction" . 

One approach to solving this problem is to derive equations of motion for q, £ and an optimal condition for 
rj and then solving a shooting problem to find n and the initial conditions for £ which allow equation ([6]) to be 
satisfied. We can derive the equations of motion by enforcing the reconstruction relation (j4|) as a constraint using 
Lagrange multipliers p q 6 T*Q. This approach leads to the following variational principle. 

Definition 4 (Variational principle for geodesic boundary value problem with symmetry). We seek 
(p, q) G T*Q and £ £ Q for t £ [0, 1], and rj £ H , which satisfy 

ss = 5 j!l { ^ A °* + ( Pq >cYt q ~ x z Gq ) T Q d< = ' (7) 

subject to 

q\t=o = q A , q\t=l=R v q B , (8) 

where we allow p q , q, £ and rj to vary. 

From this variational principle we can derive the equations of motion, which can be used in solving the shooting 
problem. Before we do this, we recall the definition of the cotangent-lifted momentum map: 

Definition 5. Given an action of a Lie algebra q on Q, the cotangent-lifted momentum map 3 : T*Q — > g is defined 
from the formula 

(J(P«r)» C> fl = (P gj -^C3>T.Q (9) 

for all £ € g. Since we have two Lie algebra actions, we shall write Jq for the cotangent-lifted momentum map 
corresponding to the left action X G of g on Q, and Jjj for the cotangent-lifted momentum map corresponding to 
the right action X H of f) on Q. 
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Lemma 6 (Equations of motion for geodesic problem). At the optimum, the following equations are satisfied 
(weakly, for appropriate pairings): 



_d 

d7 



^(T q (X^q))*p q 
M - 3 G (Pg) 



0, 

0. 

0. 



Furthermore, 



J Jj(Pg)l 



0. 



(10) 

(11) 

(12) 
(13) 



Remark 7. The end-point condition (|13j) at time t=l arises from minimising over r\ and ensures we have the 
shortest path over Q/H. 



Proof. The proof is a direct calculation. 



SS 



dq\ 



±(p q ,6^-S(X ( q) 



v dt/ T .„ 1 V<"~df 



(5l;,AZ-J G (p q )) a +(dp q ,^- 



T'Q 



dt 

T'Q 

^ + {T q ( X - q )Yp q M 



) dt + 


(Pq,Sq) T , Q 


T*Q 





1 

t=0 



Since Sp, Sq and <5£ are all arbitrary, stationarity 5S = implies equations (|10til2|> and their appropriate pairings. 
The boundary term becomes 



(P, Sq) 



T*Q 



1 

*=0 



= {p,T v (R n q) ■ Sr]) T . Q 
= (P,^ H q) T , Q 
= (Jff(Pg).7> 6 



where 7 is the generator of Sn, i.e. 



5n 



_d 

7k 



e=0 



cxp(e7)?7. 



Hence, we also obtain equation (|13p and its appropriate pairing. 



□ 



Lemma |5] reformulates the geodesic calculation as a shooting problem in which one seeks initial conditions for 
p q such that f/^i = R v q B where r\ is fixed by the condition (|13[) . Next we show conservation of the momentum 
map Jh', this will enable us to transfer condition (|13p from t — 1 to t — 0. 

Lemma 8 (Noether's theorem for geodesic problem). The system of equations flltilty) has a conserved 
momentum 3fj(p q ). 

Proof. The problem in Definition [?] is invariant under transformations 

q — > i? Q <Z, a E H, 

which are generated by 7 g f). This means that the variational principle in Definition U is invariant under application 
of the cotangent lift (i.e., the dual of the inverse of its infinitesimal transformation in Q), namely 

5q = X?q, 5p q = -(T q (X^q))*p q , 6w = 0, 

where for convenience we have inserted the time dependence 

7 = if i <*<l- 
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Substitution of this infinitesimal transformation into the variational principle gives 



o = ss 



dt 



{p q ,X„ H q) T , 



=0 
t=to 

t=0 



(7, Jff(Pg)> 6 



t=t 



T'Q 



dt 



p q +(T q {x£))*p q ,6q 



T*Q. 



dt 



(j> q ,Sq) T , 



n t=to 



t=o 



Since this equation holds for any < to < 1, it follows that the quantity 3n(Pq) is conserved. □ 

Combining this conservation result with equation (|13|) gives the following easy corollary. 

Corollary 9 ( Vanishing momentum). The conserved momentum satisfies J(p q ) = for all times < t < 1. 

Proof. Lemma |6] states that this quantity vanishes for t = 1, and Lemma [8] states that it is conserved; hence, it 
always vanishes. □ 

Corollary [9] implies that solutions of the optimal control problem all have vanishing right action momentum 
map Jff(p g ) = 0. This is what facilitates the "un-reduction" . Namely, we can compute on Q instead of Q/H by 
keeping Jjj(p g ) = 0. To obtain the shortest path between two points in Q/H by solving in Q, select a point q € Q 
which is a member of the equivalence class which is the initial point in Q/H, and find initial conditions for p q such 
that JniPq) = 0; so that the solution to equations (|10fll2[) satisfies q\t=i = R v q for some n e H. Computationally, 
there are reasons why solving the problem in this form may be difficult. In Section 13.31 we shall describe how the 
difficulty arises for the curve matching problem specified in the Introduction. In this paper, we shall introduce a 
reformulation of the problem for which there is a single fixed value of q\t=i- 

Before introducing the reformulation, we define the ad and ad* operations for the Lie algebra g and briefly 
discuss the reduced equation for the Lie algebra variable ^ e g. The latter is the Euler-Poincare equation for 
Hamilton's principle with Lagrangian given by the energy (^,^4^) /2, where A : g — > g* is the positive-definite 
self-adjoint operator in Definition [2] of the geodesic matching problem. 

Definition 10 (Notation for the ad and ad* operations). We define the operation ad : g x g — > g as 

— adtj 7 = [lj, 7] = uij — ju, uj , 7 e g, 



and define its dual ad* : g x g* — > g* as 

(ad* /x,7) fl = (^,ad w 7) B , /1 e g. 

Lemma 11 (Reduced equation for geodesic problem). The Lie algebra variable £ for the 
problem stated in Definition^ satisfies 

—M + ad* A£, = 0, 
dt s 

weakly, in the sense of the pairing ( ■ , ■ } g : g x g* — > R. 



esic matching 
(14) 
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Proof. For any 7 G Q, we have, upon substituting equation (|12[) . 



d , . , \ d , 

f^) T , Q + (^^- d i) T * Q 

= (p*, -T q {Xfq) ■ X°q + T q (Xy q) • Xf q) T . Q 

= {p^ x b,m) T *Q 

= -( 7 ,ad|^) fl . 

Consequently, we obtain the result stated, since 7 is an arbitrary clement of g. □ 

We will next define a modification of the problem stated in Definition [2 which has the advantage that the 
endpoint conditions do not contain a free reparameterisation variable. This reformulation is more amenable when 
solving the curve matching problem numerically, for example. We shall proceed to show that solutions of the 
modified problem can be transformed into solutions of the problem stated in Definition [2 

Definition 12 (Reparameterised geodesic problem with symmetry). Let Q be a manifold, let G be a Lie 
group acting on Q from the left, and let H be a (possibly different) Lie group acting on Q from the right that 
commutes with the left action of G, with corresponding Lie algebras Q and f) ; and corresponding Lie algebra actions 
X G and X H respectively. Furthermore, let A : g — > Q* be a positive-definite self-adjoint operator. We seek 

• a one parameter family q of points on Q for t £ [0, 1], 

• a one parameter family £ of elements of g for t G [0, 1], and 

• v G t), 

which minimise 

1 1 

z 

where (•, -} fl is the usual inner product on Q , subject to the constraints 

[Reconstruction relation] — q = X^q + X^q, (15) 

[Initial state (Template)] q t=0 — q , (16) 
[Final state (Target)] q\ t =i = q B , (17) 

and q A , q B are chosen points on Q. 

Remark 13. Note that in this modified definition, we do specify the final boundary condition for q without allowing 
arbitrary symmetry transformations using H . However, we also introduce an additional variable v which moves q 
in the direction of symmetries generated by f) . 

We shall derive the equations of motion associated with this modified problem, and the associated conservation 
laws. These will lead us to conclude that it possible to construct solutions of the problem in Definition [2] out of 
solutions of the problem in Definition 1 121 and the latter can be solved as a shooting problem in which the boundary 
conditions are explicitly specified, rather than as an algebraic condition. As before, we can derive the equations of 
motion for q, £ and the condition for v by enforcing the reconstruction relation (|15[) as a constraint using Lagrange 
multipliers p q G T^Q, leading to the following variational principle. 
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Definition 14 (Variational principle for reparameterised geodesic problem with symmetry). We seek 
(p, q) G T*Q and £ £ q for t G [0, 1], and v G f), which satisfy 

5S = S f a \ { ^ M) * + { Pq 'ld q + X t q + X " q ) T Q dt = °' (18) 

subject to 

q\t=v = q A , q\t=l=q B , (19) 

under variations of p q , q, £ anc? i/. 

Proceeding just as before, we can use variational calculus to obtain the equations of motion, as described in the 
following lemma. 

Lemma 15 (Equations of motion for reparameterised geodesic problem). At the optimum, the following 
equations are satisfied in the sense of appropriate pairings: 

-^q-Xfq-X?q = 0, (20) 

Ap g+ T q (xG q + X? q y -Pg = 0, (21) 
M-iaiPq) = 0. (22) 



Furthermore, 



Proof. 



[ J H ( Pq )dt = 0. 

Jt=o 



(23) 



SS = jf {6£, M) B + (s Pq , ^ - (Af + X? ) + (p q , 5^ - {X° + Xg) q - T q (Xfq + X v H q) 5q} dt 



SV, / J H (Pq)dt 

Jo 



h 

Since Sp, 5q and 6£ are all arbitrary we obtain equations (|20M23jl . □ 



Proceeding as before, we can transform (|23|) into an initial condition by making use of the conservation of the 
right-action momentum map, J# . 

Lemma 16 (Noether's theorem for reparameterised geodesic problem). The system of equations §21V2 6 2\) 
has a conserved momentum Jh ■ 

Proof. The problem in Definition fT?] is invariant under transformations 

q — > qct, a G H, 

which are generated by 7 =G f). This means that the variational principle in Definition fT?] is invariant under 
application of the cotangent lift (i.e., the dual of the inverse of its infinitesimal transformation in Q) namely 

5q = X?q, 5p=-(T q (X*q))*p, 5u = 0, 
where for convenience we have inserted the time dependence 

7 = if t <t<l. 
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Substitution of this infinitesimal transformation into the variational principle gives 

fto 



= SS 



>, - X^q - X?q, T*Q^J (p, S^q - (T q {xf ' q) + T q (X°q)) ■ Sq ^ 



dt 



to 



$P, ~^q - Xf q - X?q 



=0 



(7, JH(Pq)) f 



t=0 
t=t 

t=0 



=0 



T'Q. 



dt- 



(p, $q) 



T'Q 



t=to 
t=0 



Since this equation holds for any < to < 1, then the right-action momentum map Jh is conserved. 
Corollary 17 ( Vanishing momentum). The conserved momentum satisfies Jh = for all times t. 
Proof. After noting that Jh is conserved, Lemma [15] gives 



□ 







□ 



Next we show that £ obtained from Definition [14] satisfies the same Euler-Poincare equation as £ obtained from 
Definition 

Lemma 18 (Reduced equation for geodesic problem). The Lie algebra variable £ obeys equation (JT5J) weakly, 
i.e., for an appropriate pairing. 



Proof. For any 7 £ g, wc have, upon substituting equation ([22 

d 



dt 
d_ 

dt 



( Pq ,X^q) TtQ 
dPq X G n\ 



Pq ,T q (Xfq) 



dq\ 

dt/ T , Q 



- (T q (Xfq + X?q)Y ■ P q ,^l) T , Q + <P„ T q (X?q) ■ (Xfq + X? q)) 
(p q , -T q {X<?q + X?q) ■ X°q + T q (X°q) ■ (Xfq + X*q)) TtQ 



T'Q 



p q , Xg tQ q)^ n + (p g , T q (X°q) ■ X?q - T q (X?q) ■ X^q) 



T'Q 



7 ^/T'Q 



=0 



= <[7,£],Jg(p„)> 6 

= -( 7 ,ad|^) 1) , 

and we obtain the stated result, since 7 is an arbitrary clement of g. Here, the undcrbraced term vanishes since the 
left and right group actions commute with each other. □ 

This means that we can show that the two problems produce equivalent solutions provided that the initial 
conditions for £ are the same in both cases. The following theorem establishes this result. 
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Theorem 19. Let q, p^, v , £ be obtained from the solution of equations i2(M23\) , and define 

tp = cxp(vt) G H. 

Then the transformed variables constructed from 

q = R 4 ,-iq, Pq =T^(p 7[ ), £ = l ie[0,l], (24) 

^ i.e. the cotangent lift of tp) satisfy equations U(Mlty) together with the boundary conditions i(5I61 with n = ipt"2i- 
Hence, q, r] and £ form a (local) extremum for the problem in Definition^ 

Proof. First we take £ and £ from equations (|12j) and (|2"2"|) respectively, and show that £ = £. Since the left and 
right actions commute, we have that 



T*M 



M 



= {Pq,T q R, p X 1 q) Tm 
= (T* q R 4 ,{p 7i ),X^q) T , M 

= (7,^> , 

and hence £ = £. Next we verify the equations for q and p q . Taking the time derivative of q, we have 

di q = dt {R *' q) 

= T q Rip ■ —q + X" (R^q) , 



and so 



d _ n-1 ( d ^ vH^ 



d - 5 = (W-^-^ 

= (TgR^y 1 Xj~q 



= Xfq 
= Xfq, 

as required. To check the time evolution equation for p q , we take the inner product with an arbitrary tangent vector 
v G T q Q, to find: 

(T g i? v ,) • + (pg, (TqRip) ■ v ^ 

T* (xfq + X?q) ■ ^ (T,i^) • u)^ 

+ (g ? ,T f (^ r 5)-(T g ii ).t;> T . g 
P\,T q (xfq) ■(T q R i ,)-v) T ^ 

TZRM^fxfqj-v) 



T*Q 



= - T 



Pq ,T q (xfq)- V ) TrQ 



Cotter and Holm 



Geodesic boundary value problems with symmetry 



December 21, 2009 10 



as required. 

It remains to check the boundary conditions. Trivially, q\t=o = q A , q\t=i = o\t=i4 l \t=i = Q Br l: as required. 
Finally, we need to check the end condition (fT3|) . From Corollary |T7l we have 

J#(p g )|t=o = Jjj(%)|t=o = 0, 
and Lemma [8] implies that Jh {p q )\t=i = 0. Hence the boundary conditions are satisfied. □ 



3 Examples 

In this section we describe examples of the reparameterised geodesic problem with symmetry, and discuss its 
applications to the characterisation of the shape of curves and surfaces. 

3.1 Example: SO(3) 

We illustrate our results with the case of the action of 5*0(3) on itself which gives rise to the equations of a rotating 
rigid body. We consider the problem in which the end point boundary condition is only determined up to a rotation 
of the rigid body about its z-axis. Of course, this problem can also be solved by picking reduced coordinates, but 
we use it as here as a simple example. 

Definition 20 (Optimal control of a symmetric rigid body). Let Q(t) be a one-parameter family of matrices 
in 50(3). Let u(t) be a one-parameter family of matrices in so(3). Let Rg be a rotation in the z-axis by an angle 
9. We seek Q(t), uo(t), and Rg which satisfy 

■ r 11 

mm / 

",0 Jo 

subject to the constraints 



min / — ( w, Iu> ) dt 

Ulfi Jq 2 \ I SO(3)X«io(3)* 



[Reconstruction relation] Q(t) = oj(t)Q(t), (25) 
[Initial (Template)] Q(0) = Q°, (26) 
[Final (Target)] Q(l) = Q 1 R g , (27) 



where L is a chosen symmetric matrix. The dynamical constraint H25\) allows the reconstruction of the curve 
Q{t) £E SO(3) on the Lie group from the optimal right- invariant (spatial) angular frequency 

w{t) = QQ- l {t), 

in the Lie algebra so(3). The other constraints specify the starting and ending points of the curve Q(t) £ SO(3). 

This problem is an example of the optimal control problem in Definition [21 with the manifold Q being 5*0(3), 
the group G being 50(3) acting from the left, and the group H being 50(2) acting from the right. We identify 

q = Q, w = £, Rg = n, A = L and p q = P, 

where P € T*50(3) is the conjugate momentum to Q. Application of Lemma [6] gives the following dynamical 
equations: 

Q-ujQ = 0, (28) 
P + lu t P = 0, (29) 

Llu + ^(PQ t -QP T ) = 0, (30) 

which are the equations for a rotating rigid body. The last line follows from the definition of the left-action 
momentum map for 50(3), namely 

(J L (P)M.oW = (^-^0) T , SO (3) = -tr(P T SuQ) = -tr((PQ T ) T ^) = (PQ T ,~5u;) so{3) , 
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where 5u> is an arbitrary antisymmetric 3x3 matrix. (Hence, the antisymmetric combination in equation (|30p .) 
The end point condition (which comes from minimising over 9) becomes 

(P,Qw) T , so(3) \ t=1 =0, (31) 

where 

/ 1 0\ 

(32) 

and Lemma |9] implies that the quantity (P,Qw) T , SO pj (which is the z-component of the angular momentum) is 
zero for all times t. Lemma [14] states that uj satisfies the Eulcr equations for a rigid body: 

^(Iw) + = 0, 

where we define 

ad^7 = [w,7] = cc>7 - ju, cj,7<Gso(3), 

and 

(ad* M,7> = (^,ad w 7) , p £ so(3)* . 

For this problem, obtaining a solution is simple, since one can define coordinates on T(50(3)), and remove 
the coordinates associated with the Rg direction and the corresponding vanishing conserved momentum, and solve 
a two-part boundary problem for the remaining coordinates. However, we wish to develop a methodology for 
numerical discretisations of infinite-dimensional problems where it is less clear how to do this. Hence, we define the 
following formulation which makes use of a time-varying "relabelling" transformation in the Rg direction. Theorem 
[19] states that to obtain solutions to equations (|28H30|) . we can solve the following modified problem: 

Definition 21 (Reparameterised optimal control of a symmetric rigid body). Let Q(t) be a one-parameter 
family of matrices in 50(3). Let uj(t) be a one-parameter family of matrices in so(3). Let v be the generator of a 
rotation about the z-axis, which may be written in the form 

v = 6w, 



where w is defined in equation (|32j) . and 9 £ 
We seek Q(t), U(t) 7 and v which satisfy 



mm 



— ((jj.Iu) dt 



ui-M Jq 2 \ / so(3)xso(3)* 

subject to the constraints 

[Reconstruction relation] Q(t) = uJ(t)Q(t) + ~Q{t)v, (33) 

[Initial (Template)] Q(0) = Q°, (34) 

[Final (Target)] Q(l) = Q\ (35) 

where I is a chosen symmetric matrix. 

Lemma [15] states that the solution to this problem satisfies the following equations: 

ti-uQ + Qv = 0, (36) 
T + UJ T P -T>v T = 0, (37) 

IuJ+^(PQ T -QP T ) = 0, (38) 

with end-point condition 

<^Q^> T .so(3)l*=i = °- ( 39 ) 

Lemma [T71 states that 

(P,Qw) T , so{3) =0, (40) 
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for all t. 

Hence, to obtain a solution to equations (|28H31[) . we solve the two-point boundary value problem given by 
equations (|36H37[) with boundary conditions (|34II35[) . This can be formulated as a shooting problem, in which we 
seek v (or, equivalently, 8) and initial conditions for P satisfying equation (|40|). such that the end point boundary 
condition (|35|) is satisfied. We then construct the reparameterisation matrix R(t) from 

R(t) = exp(>t), 

and use equation (|24p to reconstruct the solution in the form: 

Q(t) =Q{t)R(t) and P(t) = Q{t)R T (t), implying u(t)=W(t), 

since, e.g., R T = R^ 1 implies PQ = PQ T . Then substituting these relations into equations (|28ti30|) and equation 
([31]) recovers equations (|36M38j) and equation (j40|) . 



3.2 Example: SE(3) 

We next describe the example of the action of SE(3) on itself from the left, with 5*0(2) acting from the right. 
This example could describe a docking problem of a spacecraft onto a space station. The spacecraft can apply 
torque to rotate about a central point, or can apply thrust to move itself in the direction in which it is pointing, 
and we wish to dock the spacecraft using minimal energy. In the language of image registration, this is known as 
rigid registration. We consider the problem in which the end point boundary condition is only determined up to 
a rotation of the rigid body about its z-axis. In the spacecraft analogy, this corresponds to a docking procedure 
which does not require the spacecraft to have any particular orientation about the z-axis when docking. As in the 
previous example, this problem can also be solved by picking reduced coordinates, but it serves as a prototype for 
infinite dimensional problems. 

Following the notation of [Hol09], we represent an element q of SE(3) as a 4 x 4 matrix: 

Q r 
1 

where Q is an orthogonal matrix, and r £ R 3 . We represent an element p q of T*SE(3) as a 4 x 4 matrix: 



Pq 



_ (P P 

o or 



where P G T*SO(3) and p £ M 3 . Finally, we represent an element of the corresponding Lie algebra se(3) as another 
4x4 matrix: 

<-(Z o). 

where u> is an antisymmetric matrix, and o£l 3 . The reconstruction relation is then given by 

1 = X t 1 = & = ( o 



We write the energy cost for the system as 



1 



S= Edt=- (Z,M) st(3) dt, 
Jt=o z Jt=o 



where A is a symmetric positive definite 4x4 matrix given by 
The starting condition is specified as 



I b 

b T c 



q\t=o 



Q A r AN 
1 I ' 
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which specifies the starting orientation and position, and the end condition is specified as 

(Q B Re r E 

«l*=i= o i 

where Re is a rotation through any angle 9 about the z-axis. If we solve the problem in Definition [2l then Lemma 
[6] gives the dynamical equations 



Q 


= uQ, 


(41) 


r 


= ojr + v, 


(42) 


P 


= -^ T P, 


(43) 


P 


T 

= u p, 


(44) 




= o, 


(45) 


P 


= o, 


(46) 



SE_ 

where 8E/Suj = Iu> + bv T — vb T , SE/Sv = cv + 2u>h. The corresponding end point condition is 

(P,Qw) T . SO (3)lt=i=0, (47) 

where w is defined in equation (f3"2"|) . as for the SO(3) case. From LemmalU this quantity vanishes for all t. Theorem 
[Tnithen states that a solution to these equations can be obtained by solving the following reparameterised equations: 

~Q = ZJQ + aQw, (48) 
r = ZJr + v, (49) 

P = -uF~P-aQw T , (50) 
P = u T p, (51) 

S E 1 T T 

-8=+2 {PQ QP) = °' (52) 

5E , . 

t^+P = 0, 53 
ov 

with a£i Q\t=a = Q A : r\t=o = rA '■ Q\t=i = Q B , r\t=i = r B , an d endpoint condition 

(P,Otw) T . sf O(3)l*=o=0- ( 54 ) 

This gives a two-point boundary value problem with a constraint on the initial conditions plus an extra parameter, 
which can be solved as a shooting problem by finding a and P (subject to the constraint) such that Q and r reach 
their target values Q B and r B . A solution to the problem in Definition [2] can then be reconstructed by defining 
R{t) = exp(awt), and using the following formulae: 

Q(t) = Q(t)R(t), P(t)=P(t)R T (t), r(t)=r(t), p(t) =p(t), w(t)=57(t), v(t) =v(t). 

Substituting these relations into equations (|41M6p and equation (|47|) recovers equations (|48H53p and equation |54|) . 

3.3 Curve matching 

In this section we return to the problem described in Definition [TJ and discuss a number of practical issues which are 
addressed by the formulation discussed in this paper. The aim of solving the problem is to find a characterisation 
of the simple closed curve T B in terms of the reference simple closed curve T A , together with a scalar periodic 
function p(s) which specifies the initial conditions for the normal component of the conjugate momentum p(s;t). 
In this case, Q is the space Emb(5 1 ,IR 2 ) of functions q : S 1 — > R 2 , G is the group Diff(M 2 ) of diffeomorphisms of 
]R 2 which acts on Q from the left 

and H is the group Diff(S' 1 ) of diffeomorphisms of S 1 , which acts on Q from the right 

$>R(g,qKs) = q(g(s)), VsES 1 . 
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The left and right actions can be expressed succinctly as 

GQ = Emb(S\ G • R 2 ) , HQ = Emb(# ■ S\R 2 ). 

It is clear from this that the actions of G and H commute with each other. 
Lemma [5] gives the dynamical equations 

d 

— q(s;t) = u{q(s;t),t), 

^P(s;t) = -(Vu(q(s;t),t)) T ■ p{s;t) , 
where the velocity u is defined weakly from the equation 

(w,Au(-,t)) x = / p(s;t) ■ w(q(s;t),t)ds, 
Js 1 

where (-, -) x is the inner product on the vector fields X(fi) associated with the norm || ■ \\x, for any test function 
w E X* . This equation has the weak solution 

m(x,t) = ^2p(t,s)S(x - q(t,s)) , (55) 

which is the singular-solution momentum map, Jsing discussed in |HM04| . The end condition is 

p(a; 1) ~q(«; 1) = 0, Vs £ S\ (56) 

and Lemma [8] states that this conserved momentum vanishes for all t. This corresponds to p(s; t) being normal 
to the curve parameterised by q(s;t). Hence, to find geodesies between T A and T B , we solve a shooting problem 
and seek initial conditions p(s;0) with vanishing tangential component, which fix q(s; 1) = q B (n(s)) for some 
n £ Diff (S 1 ). Having solved this problem, one can describe T B entirely in terms of p(s) = p(s; 0) • n(s) where n is 
the normal to T A . The solution to the problem also provides the distance between the two curves. 

There are a number of difficulties with solving such a shooting problem numerically. The parameterisation of 
the curve must necessarily be discretised numerically, typically by a list of points, as in MM06, Cot08j, which can 
be obtained from a piecewise-constant representation of q as a function of s |Via09j , or as piecewise linear geometric 
currents |VG05| . Having taken the discretisation, the reparametcrisation symmetry is broken (although a remnant 
of it is left behind, as described in |Cot08| ) which means that it is difficult to obtain a discrete form of the end 
condition q(s; 1) = q B (i](s)). As described in the Introduction, this problem has been approached by proposing 
various functionals which arc minimised when q(s; 1) and q(s) overlap the most. However, these functionals produce 
quite a complicated landscape with local minima, and the case of studying large deformations we have found that 
they can result in odd artifacts in the shooting process. Also, the changes in these functionals with respect to 
measurement error are quite technical to quantify which makes statistical inference more complicated. 



Another difficulty is that of adaptivity. As illustrated in Figure [TJ constraining p to be normal to the curve 
means that any local large deformations give rise to large amounts of stretching which then results in loss of accuracy 
in the approximation of the functional used to enforce the end condition for q. One possible way to avoid this is to 
adaptivcly refine the grid point density in the initial curve during the shooting process. 

These difficulties are all removed if, instead, one solves the following problem: 

Definition 22 (Reparameterised curve matching problem). Letq(s;t) be a one-parameter family of parame- 
terised curves in M 2 , with s £ [0, 1] being the curve parameter and t £ [0, 1] being the parameter for the family. Let 
u(x;t) be a one-parameter family of vector fields on M 2 . Let v be a vector field on S . We seekq, u and v which 
satisfy 




subject to the constraints 

d _ d 

[Reconstruction relation] —q(s;t) = u(q(s;t),t) + v(s)—q(s;t), (57) 

Initial state (Template)] q(s;0) = q A (s), (58) 
[Final state (Target)] q{s; 1) = g B (s), (59) 

(60) 
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Figure 1: Figure illustrating the way in which deformation takes place in the parametcrisation- independent geodesic 
equations for embedded curves. The initial curve is shown on the left, and the final curve is shown on the right. 
Since the momentum is constrained to be normal to the curve, and since the velocity is obtained by applying a 
smoothing kernel to the momentum, the change in the shape emerges as local stretching of the curve, and the 
discrete points defining the shape become separated. 



where \\ ■ \\v is the chosen norm which defines the space of vector fields V. 

Lemmas [15] and [17] state that extrema of this problem can be obtained by solving the shooting problem 

d _ d _ 

— q{s;t) = u(q(s;t),t) + v(s)—q(s;t), 
at os 



with boundary conditions 



j-p{s;t) = -p(s;t)-Vu(q{s;t),t)- — {v(s)p{s;t)), 
Ot os 

{w,u(-,t)) v = / p(s;t) ■ w(q(s;t),t)ds, \/w G V* , 
Js 1 



q(s;0) = q A (s), p(s; 0) • |V = 0, q(s; 1) = q B (s). 

os 



The aim is to find v(s) and normal components of p(s; 0) such that these boundary conditions are satisfied. Note that 
in this modified problem, the boundary condition for q(s; 1) is specified pointwisc (i.e. there is no reparamcterisation 
variable i] in the boundary condition). This means that the error can be described directly in terms of 



\\q(-,l)-q B \\ 2 

for some chosen norm, which can be discussed in terms of measurement errors directly. 

Theorem [19] then states that a solution to the problem described in Definition [T] can be reconstructed via 

d 

q(s;t) = q(n(s;t),t), p{s;t) = p(n(s;t),t)—n(s;t) 

where 

r)(s;t) = exp(i/(s)£). 

This transformation produces a equivalent shooting problem in which the end condition for q is now fixed. 



4 Summary and outlook 

In this paper, we studied an optimal control problem on a Lie group in which the end boundary condition is 
specified only up to a symmetry. We showed how this problem can be transformed into a modified problem in 
which the end boundary condition is fixed, but an extra parameter is introduced in the reconstruction relation, 
and proved that the two problems are equivalent. This approach is motivated by the problem of computing a 
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geodesic on the diffcomorphism group in the plane which takes one curve to another. The transformation gives 
rise to a system of equations for a parameterised curve in which the end boundary condition for each value of 
the parameter is fixed. This means that when a discrete approximation of the curve is used to solve this problem 
numerically, the end boundary condition can again be fixed exactly. In particular, when solving a shooting problem, 
this means that the error between the computed curve and the target curve can be computed simply by measuring 
the Euclidean distance between points for each value. This method extends straightforwardly to the problem of 
finding geodesies in the three-dimensional diffcomorphism group which take one parameterised surface to another, 
with the end boundary condition specified only up to reparameterisations of the target surface. This problem has 
many applications in, for example, biomedicine, since it allows topologically equivalent surfaces to be characterised 
by a momentum field distributed on the template surface. Such momentum fields exist in a linear space and so can 
be manipulated using linear techniques and still a topologically equivalent surface will always be obtained. 

In the standard approach to planar image registration, the problem of registering a specified closed curve 
(called the template) at time t = to another (the target) at time t = 1 amounts to deforming the space M 2 in 
which the template curve is embedded until it matches the fixed target image to within a certain tolerance. Here, 
we considered a manifold Q comprising the space of closed curves S 1 embedded into the plane IR 2 , written as 
Q = Emb(S' 1 , IR 2 ). There are two Lie group actions available for manipulating the closed curves in this description. 
The action G x IR 2 — > IR 2 of the Lie group G = Diff(IR 2 ) by composition from the left deforms the range space 
M 2 , and thereby drags along a curve embedded in it as GQ = Emb(5 1 ,G • R 2 ). This left action produced the 
singular-solution momentum map, Jsing in equation (|55[) . which introduces the paramcterisation of the closed curve 
by its position and momentum supported on a delta function defining the curve in IR 2 . Under this left action of 
G, the curve preserves the initial parameterisation of its domain space in S , although the current positions of the 
S 1 labels in the plane IR 2 will change as the range space is transformed. Alternatively, the action H x S 1 — > S 1 
of the Lie group H = Diff(S' 1 ) by composition from the right transforms coordinates in the domain space 5" 1 as 
HQ = Emb(ff • S 1 , IR 2 ), while keeping the curve fixed in the range space IR 2 . 

The present paper discussed how the left action of Diff (IR 2 ) and the right action of Diff (S 1 ) on Q = Emb(S' 1 , IR 2 ) 
may complement each other in formulating a variational approach for registration of curves under large deformations. 
The left action of Diff (IR 2 ) corresponds to deforming the curve by a time-dependent transformation of the coordinate 
system in which it is embedded, while leaving the parameterisation of the curve invariant. The dynamics of this 
deformation is formulated as an Euler-Poincare equation for Jsing G X(IR 2 )* that results in Hamilton's canonical 
equations for the momentum and position variables of the curve that comprise the singular-solution momentum map 
(1551) . This solution provides the dynamics for curves that fulfills D'Arcy Thompson's vision of shape transformation 
|Thol7j . This vision underlies common practice in image registration |Beg03| . 

The right action of Diff(S' 1 ) corresponds to adaptively reparameterising the S 1 domain coordinates of the 
curve. This rcparametcrisation of the curve could be useful, for example, in designing numerical methods that 
enhance the resolution of its main features as it deforms in M 2 . As we have seen, this right action unlocks the 
parameterisation in the control problem to allow it more freedom for matching the curve shapes using an equivalent 
boundary value problem, without being constrained to match corresponding points along the template and target 
curves at the endpoint in time. As explained above, the action of Diff(5 1 ) from the right gives us the momentum 
map Js : T* Emb(S' 1 ,IR 2 ) — * X(S)* , which we used to ensure that the momentum of the curve has no tangential 
component. This momentum map for right action is given explicitly as 



os 



The two momentum maps may be assembled into a single figure as in |HM04| : 



T* Emb(5 





X(R 2 )* 



X(S 1 )* 



We use the compatibility of these two momentum maps proven in |HM04j to divide the curve matching problem 
into independent registration and reparameterisation problems, leading to the reformulation of the curve matching 
problem as an equivalent geodesic boundary value problem. 
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We are currently developing numerical algorithms for the transformed equations applied to embedded curves 
and surfaces. As noted in |Via09j . applying a piecewise constant representation to the q and p variables in the 
untransformed equations results in a set of ordered points. When this approach is extended to the transformed 
equations, a finite volume method is obtained, with the extra terms taking the form of an advection term in the q 
equation, and a continuity term in the p equation, which are very well developed in the finite volume approach. We 
will also investigate discontinuous higher-order polynomial representations of p and q which lead to discontinous 
Galcrkin methods. Since the error in the transformed problem can be quantified in terms of the Euclidean distance 
between points on the curve for each parameter value, the reparameteriscd formulation also makes it much easier 
to perform Bayesian studies in which one observes points on a curve with observation error from some probability 
distribution, and then one attempts to estimate the probability distribution for the initial conditions of p (and v) 
for which specified points on the curve match the actual position of the observed points, after applying the time-1 
flow map of the transformed geodesic equations for p and q. This approach provides a considerably simplified 
observation operator to which algorithms such as the Monte Carlo Markov Chain algorithm can be applied. 
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